package com.sr.srinterpolation.util;

import android.util.Log;

/**
 * 牛顿插值法
 */
public class NewtonInterpolationUtil {
    private static final String TAG = "NewtonInterpolationUtil";

    /**
     * 返回牛顿插值多项式
     * 牛顿插值多项式：y = c0 + c1(x - x0) + c2(x - x0)(x - x1) + ...
     *
     * @param x {1, 2, 3, 4, 5, ... , n}
     * @param y {1, 2, 4, 8, 16, ... , 114514}
     * @return {{分子数组}{分母数组}}
     */
    public static int[][] polynomial(int[] x, int[] y) {
        int n = x.length;
        // cn的分子和分母数组
        int[] cNums = new int[n];
        int[] cDens = new int[n];
        // 遍历求cn
        cNums[0] = y[0];
        cDens[0] = 1;
        cDens[1] = 1;
        for (int i = 1; i < n; i++) {
            int num = y[i];
            // 每一项的分母通分后的数字
            int cReduceDen;
            // 每一项的通分倍率
            int[] reduceRadios;
            // 计算以上两项
            if (i < 2) {
                cReduceDen = 1;
                reduceRadios = new int[]{1};
            } else {
                cReduceDen = Math.abs(lcm(cDens, i));
                reduceRadios = reduceRadio(cDens, cReduceDen, i);
            }
            // 分子
            num *= cReduceDen;
            // 遍历计算ci之前的每一项，yi - c0 - c1(x - x0) - c2(x - x0)(x - x1) - ... - ci-1(x - x0)...(x - xi-2)
            // xp为(x - x0)(x - x1)...(x - xj-1)部分的值
            int xp = 1;
            for (int j = 0; j < i; j++, xp *= x[i] - x[j - 1]) {
                num -= cNums[j] * xp * reduceRadios[j];
            }
            // 计算分母
            int den = cReduceDen * xp;
            // 约分并记录约分后的ci
            int gcd = MathUtil.gcd(num, den);
            cNums[i] = num / gcd;
            cDens[i] = den / gcd;
        }
        // 计算每一项的分母通分后的数字
        int cReduceDen = Math.abs(lcm(cDens, n));
        // 计算每一项的通分倍率
        int[] reduceRadios = reduceRadio(cDens, cReduceDen, n);
        // 计算分子部分
        int[] nums = nums(x, cNums, reduceRadios);
        // 约分得出最终结果
        int[][] result = MathUtil.reductionOfAFraction(nums, cReduceDen);

        for (int num : nums) {
            Log.i(TAG, "num: " + num);
        }
        for (int resultNum : result[0]) {
            Log.i(TAG, "resultNum: " + resultNum);
        }
        for (int resultDen : result[1]) {
            Log.i(TAG, "resultDen: " + resultDen);
        }
        return result;
    }

    // 计算分子部分
    private static int[] nums(int[] x, int[] cNums, int[] reduceRadios) {
        int n = x.length;
        int[] nums = new int[n];
        // (x - x0)(x - x1)...(x - xn)部分的多项式的系数
        int[] xNums = new int[n];
        xNums[0] = 1;
        nums[0] = cNums[0] * reduceRadios[0];
        // 遍历计算每一项的分子部分，y = c0 + c1(x - x0) + c2(x - x0)(x - x1) + ...
        for (int i = 1; i < n; i++) {
            // 计算(x - x0)(x - x1)...(x - xi-1)部分的多项式的系数
            xNums[i] = 1;
            for (int j = i - 1; j > 0; j--) {
                xNums[j] = xNums[j - 1] - xNums[j] * x[i - 1];
            }
            xNums[0] *= -x[i - 1];
            // 计算nums
            for (int j = 0; j <= i; j++) {
                nums[j] += xNums[j] * cNums[i] * reduceRadios[i];
            }
        }
        return nums;
    }

    // 获取通分倍率（所有分数通分后，每项所需乘的倍数）
    private static int[] reduceRadio(int[] dens, int reduceDen, int n) {
        int[] radios = new int[n];
        for (int i = 0; i < n; i++) {
            radios[i] = reduceDen / dens[i];
        }
        return radios;
    }

    // 多个数的最小公倍数
    private static int lcm(int[] a, int n) {
        int lcm = MathUtil.lcm(a[0], a[1]);
        for (int i = 2; i < n; i++) {
            lcm = MathUtil.lcm(lcm, a[i]);
        }
        return lcm;
    }
}
